# Dickstein, Ho, and Mark (2023)
# This script includes functions that are used to generate
# counterfactuals and to report counterfactual comparisons

# * # * # * # * # * # * # * # * #
# Run Counterfactual Functions #
# * # * # * # * # * # * # * # * #

RunCounterfactual <- function(
  ExplDat,
  PlanDat,
  PremVec.Init = NULL, # Initial Value of Premium vector. (for each constructed plan ID)
  behaviorShare = .01, # percent of the market that we add who are behavioral types.
  maxR = 100, #Maximum number of iterations
  tol = 10^(-3), # tolerance
  costwinsorization = 89.37, # winsorizing cost value
  pricelimit = 1000, # limit on prices
  subsidytypeSG = "Ratio", # Other options are "Fixed" or "CSCAdjust"
  subsidytypeInd = "Ratio", # Other options are "Fixed" or "CSCAdjust"
  # Ratio: subsidy is fixed % taken off.
  # Fixed: subsidy is a fixed $ taken off.
  # CSCAdjust: subsidy changes given the second cheapest silver plan.
  fixed_markup = 0
  ){ 
  
  # Pre-algorithm cleaning
  ## Ensure that all plans exist in exploded data and all options exist in plan data
  nr <- nrow(ExplDat)
  nhh <- length(unique(ExplDat$hhid))
  ExplDat <- ExplDat[(plan_id %in% unique(PlanDat$plan_id))]
  ExplDat[, n := .N, by = c("hhid", "year")]
  ExplDat <- ExplDat[n > 1]
  ExplDat$n <- NULL
  print(paste0(nr - nrow(ExplDat), " observations removed due to choices not in ind. market choice set."))
  print(paste0(nhh - length(unique(ExplDat$hhid)), " households removed due to choices not in ind. market choice set."))
  nr <- nrow(PlanDat)
  PlanDat <- PlanDat[(plan_id %in% unique(ExplDat$plan_id) | plan_id == "Outside_Option")]
  print(paste0(nr - nrow(PlanDat), " plans removed due to never existing in choice set."))
  
  ## Generate PremVec.Init if it is not given
  if(is.null(PremVec.Init)){
    PremVec.Init = ifelse(PlanDat$plan_id == "Outside_Option", 0, 1)
  }
  
  # Preliminaries before while loop starts
  ## Plug in initial price vector
  PlanDat$price <- PremVec.Init
  avepercentdiff <- 1
  r <- 0
  
  while(r < maxR & avepercentdiff > tol){
    r <- r + 1
    print(paste0("Starting step", r))
    # Step Zero: Merge price vector into exploded dataset and create individual prices:
    ExplDat <- merge(
      ExplDat,
      PlanDat[, .(plan_id, price)],
      by = "plan_id",
      all.x = T)
                     
    AnyIndMkt <- sum(ExplDat$hhid > 0)  
    AnySGMkt <- sum(ExplDat$hhid < 0)
    
    ExplDat_List <- list()
    if(AnyIndMkt) ExplDat_List$IM <- ExplDat[hhid > 0]
    if(AnySGMkt) ExplDat_List$SG <- ExplDat[hhid < 0]
    
    subsidytype_vec <- c(subsidytypeInd, subsidytypeSG)
 
    for(i in seq_along(ExplDat_List)){              
      if(subsidytype_vec[i] == "Ratio"){
        ExplDat_List[[i]][, price_hh := ratingfactor * (1 - subs) * price, ]
      } else if(subsidytype_vec[i] == "Fixed"){
        ExplDat_List[[i]][, price_hh := ifelse(hhid < 0, 
          pmax(ratingfactor * (1 - (best_guess_sumrate / 100)) * price - subs, 0), 
          pmax(ratingfactor * price - subs, 0)), ] 
        # Note: subs in this case is (1-t)*v * base rate before
      } else if(subsidytype_vec[i] == "SCSAdjust"){
        if(r == 1){    
          if(sum(ExplDat_List[[i]]$scs) != nrow(unique(ExplDat_List[[i]][, .(hhid, year), ]))){
            stop("Each household needs exactly one second cheapest silver plan")
          }
        }  
        ExplDat_List[[i]][, price_scsr := sum(scs * ratingfactor * price), by = c("hhid", "year")]     
        ExplDat_List[[i]][, subs := pmax(subs_scs0 + (price_scsr - price_scs0), 0), ] 
        ExplDat_List[[i]][, price_scsr := NULL, ]
        ExplDat_List[[i]][, price_hh := ifelse(hhid < 0, 
          pmax(ratingfactor * (1 - (best_guess_sumrate / 100)) * price - subs, 0), 
          pmax(ratingfactor * price - subs, 0)), ] 
      } else {
        stop("subsidytype must be one of: Ratio, Fixed, SCSAdjust")
      }
    }
    ExplDat <- do.call("rbind", ExplDat_List)
    
    # Step One: Compute household expected utility for each plan:
    ExplDat[, V := x_av + .5 * omega * (x_av ^ 2) - (alpha - psi) * price_hh + fixedeffect_beta0, ]
    
    # Step Two: Compute predicted plan shares (within each household)
    ExplDat[, sumexpV := sum(exp(V)), by = c("hhid", "year")]
    ExplDat[, share := exp(V) / sumexpV, ]
    
    # Step Three: Compute expected costs for household i under plan j:
    ExplDat[, costbar := pmin((x_av + omega * (x_av ^ 2)) / alpha, costwinsorization)]

    # Step Four: Compute total cost to the insurer of the plan and the new price
    CostPreAdjustment <- ExplDat[plan_id != "Outside_Option", .(
      cost = sum(share * kappa * costbar),
      ratingfactorsharesum = sum(share * ratingfactor)
    ), by = "plan_id"]
    
    PlanDat <- merge(
      PlanDat,
      CostPreAdjustment,
      by = "plan_id",
      all.x = T
    )
    
    PlanDat[, pricenew := (1 + fixed_markup) * (beta_4 * cost + fe_beta_5) / ratingfactorsharesum * (1 + behaviorShare), ]
    if(sum(is.na(PlanDat$pricenew)) != 1){ stop("Missing new price")}
    PlanDat[, pricenew := ifelse(is.na(pricenew), 0, pricenew), ]
    
    # Step Five: Test for convergence
    avepercentdiff <- mean(abs((PlanDat[pricenew != 0]$pricenew - PlanDat[pricenew != 0]$price) / PlanDat[pricenew != 0]$price))
    print(paste0("Average Percent Change is ", avepercentdiff))
    
    # Step Six: Redefine price and clean
    PlanDat <- PlanDat[, .(
      plan_id,
      beta_4,
      fe_beta_5,
      price = pmin(pricenew, pricelimit)
    )]
    
    # Cleaning
    if(r < maxR & avepercentdiff > tol){
      vars_to_keep <- c("orig_subs_id", "year", "hhid", "plan_id", "best_guess_AV_old", 
        "best_guess_netprem_old", "best_guess_sumrate", "grossprem_old", "ratingfactor", 
        "subs", "x_av", "alpha", "omega", "psi", "fixedeffect_beta0", "kappa", "sum_concurrent_risk")
        # "old" variables are only for reporting
      if(any(subsidytype_vec %in% c("SCSAdjust"))){
        vars_to_keep <- c(vars_to_keep, "scs", "price_scs0", "subs_scs0") 
      }
      
      ExplDat <- ExplDat[, vars_to_keep , with = F]
    }
  }
  
  return(CounterfactualList = list(
    ExplodedData = ExplDat,
    PlanData = PlanData
  ))
}

# * # * # * # * # * # * # * # * #
#       Table Generation        #
# * # * # * # * # * # * # * # * #

cmax <- function(x, y) {ifelse(x>y, x, y)}

generateTableOne <- function(CounterFactualResultsPre,
  CounterFactualResultsPost,
  costwinsorization = 89.37,
  SGpremsubs = F, # Are SG households given g'ment subsidies?
  subsidytypeSG = "Ratio", # Other options are "Fixed" or "CSCAdjust"
  subsidytypeInd = "Ratio" # Other options are "Fixed" or "CSCAdjust"
){
  # RESULTS TABLE ONE SKELETON
  print("Making skeleton...")
  ntablerows <- 19
  Table1 <- data.frame(
    IndPre = numeric(ntablerows),
    IndPost = numeric(ntablerows),
    IndChange = numeric(ntablerows),
    SGPre = numeric(ntablerows),
    SGPreSamplePost = numeric(ntablerows),
    SGChange = numeric(ntablerows),
    SGPost = numeric(ntablerows)
  )
  row.names(Table1) <- c(
    "Consumer Surplus ($\\alpha > \\psi$)",
    "P($\\alpha > \\psi$)",
    paste0("Cons. Sur. ($\\alpha - \\psi >", alphampsicutoff, "$)"),
    paste0("P($\\alpha - \\psi >", alphampsicutoff, "$)"),
    "$E[\\lambda_i]$",
    "Moral Hazard Spending",
    "Share Uninsured",
    "Share Bronze",
    "Share Silver",
    "Share Gold",
    "Ave (Base) Bronze Prem",
    "Ave (Base) Silver Prem",
    "Ave (Base) Gold Prem",
    "Ave (Net) Bronze Prem",
    "Ave (Net) Silver Prem",
    "Ave (Net) Gold Prem",
    "Government Spending",
    "Employer Spending",
    "N")
  
  # Creating measures for each person:
  ## Pre-Market
  print("Collapsing to household (pre)...")
  pretab <- CounterFactualResultsPre$ExplodedData[, .(
    RA = mode1(substr(plan_id, 3, 3)),
    CS = mode1(1 / (alpha - psi)) * log(sum(exp(V))),
    p_1 = ifelse(mode1(alpha - psi) > 0, 1, 0),
    p_2 = ifelse(mode1(alpha - psi) > alphampsicutoff, 1, 0),
    PUninsured = sum(share * as.numeric(x_av == 0)), 
    PBronze = sum(share * as.numeric(x_av == 0.6)), 
    PSilver = sum(share * as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94))),
    PGold = sum(share * as.numeric(x_av == 0.8)),
    ExpAV = sum(share * x_av), 
    ExpLambda = sum(share * pmin(1 / alpha, costwinsorization)), 
    ExpMHS = sum(share * pmin((x_av + omega * (x_av ^ 2)) / alpha, costwinsorization)) - sum(share * pmin(x_av / alpha, costwinsorization)), 
    AveBronzePremNet = sum(as.numeric(x_av == 0.6) * ifelse(is.na(price_hh / ratingfactor), 0, price_hh / ratingfactor)) / sum(as.numeric(x_av == 0.6)), 
    AveSilverPremNet = sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94)) * ifelse(is.na(price_hh / ratingfactor), 0, price_hh / ratingfactor)) / sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94))), 
    AveGoldPremNet = sum(as.numeric(x_av == 0.8) * ifelse(is.na(price_hh / ratingfactor), 0, price_hh / ratingfactor)) / sum(as.numeric(x_av == 0.8)), 
    AveBronzeBaseGross = sum(as.numeric(x_av == 0.6) * ifelse(is.na(price), 0, price)) / sum(as.numeric(x_av == 0.6)),
    AveSilverBaseGross = sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94)) * ifelse(is.na(price), 0, price)) / sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94))), 
    AveGoldBaseGross = sum(as.numeric(x_av == 0.8) * ifelse(is.na(price), 0, price)) / sum(as.numeric(x_av == 0.8)), 
    GovmntSpend = ifelse(hhid > 0, 
      sum(share * subs * ratingfactor * price),
      sum(share * (best_guess_sumrate / 100) * ratingfactor * price)), 
    EmployerSpend = sum(share *  (subs / ( 1 - subs)) * price_hh)
  ),
  by = c("hhid", "year")]
  
  pretab[, Market := ifelse(hhid < 0, "SmallGroup", "Individual")]
  pretab[, CS_1 := ifelse(p_1 == 1, CS, NA), ]
  pretab[, CS_2 := ifelse(p_2 == 1, CS, NA), ]
  
  ## Post-Market
  print("Collapsing to household (post)...")
  posttab <- CounterFactualResultsPost$ExplodedData[, .(
    RA = mode1(substr(plan_id, 3, 3)),
    CS = mode1(1 / (alpha - psi)) * log(sum(exp(V))),
    p_1 = ifelse(mode1(alpha - psi) > 0, 1, 0),
    p_2 = ifelse(mode1(alpha - psi) > alphampsicutoff, 1, 0),
    PUninsured = sum(share * as.numeric(x_av == 0)), 
    PBronze = sum(share * as.numeric(x_av == 0.6)), 
    PSilver = sum(share * as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94))),
    PGold = sum(share * as.numeric(x_av == 0.8)),
    ExpAV = sum(share * x_av), 
    ExpLambda = sum(share * pmin(1 / alpha, costwinsorization)), 
    ExpMHS = sum(share * pmin((x_av + omega * (x_av ^ 2)) / alpha, costwinsorization)) - sum(share * pmin(x_av / alpha, costwinsorization)), 
    AveBronzePremNet = sum(as.numeric(x_av == 0.6) * ifelse(is.na(price_hh / ratingfactor), 0, price_hh / ratingfactor)) / sum(as.numeric(x_av == 0.6)), 
    AveSilverPremNet = sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94)) * ifelse(is.na(price_hh / ratingfactor), 0, price_hh / ratingfactor)) / sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94))), 
    AveGoldPremNet = sum(as.numeric(x_av == 0.8) * ifelse(is.na(price_hh / ratingfactor), 0, price_hh / ratingfactor)) / sum(as.numeric(x_av == 0.8)), 
    AveBronzeBaseGross = sum(as.numeric(x_av == 0.6) * ifelse(is.na(price), 0, price)) / sum(as.numeric(x_av == 0.6)),
    AveSilverBaseGross = sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94)) * ifelse(is.na(price), 0, price)) / sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94))), 
    AveGoldBaseGross = sum(as.numeric(x_av == 0.8) * ifelse(is.na(price), 0, price)) / sum(as.numeric(x_av == 0.8)), 
    GovmntSpend = ifelse(SGpremsubs == T | hhid > 0, 
      ifelse(subsidytypeInd == "Ratio",
        sum(share * subs * ratingfactor * price), 
        sum(share * subs)),
      sum(share * cmax((best_guess_sumrate / 100) * ratingfactor * price , nu * (best_guess_sumrate  / 100) * (grossprem_old / 100)))),
    EmployerSpend = ifelse(SGpremsubs == T | hhid > 0, 
      0,
      ifelse(subsidytypeSG == "Ratio", 
        sum(share * (nu / ( 1 - nu)) * price_hh), 
        subs)),
    InPreSample = mode1(!is.na(grossprem_old))
  ),
  by = c("hhid", "year")]
  
  posttab[, Market := ifelse(hhid < 0, "SmallGroup", "Individual")]
  posttab[, CS_1 := ifelse(p_1 == 1, CS, NA), ]
  posttab[, CS_2 := ifelse(p_2 == 1, CS, NA), ]
  
  ## Small Group Pre (observables)
  print("Collapsing to household (Small Group Pre)...")
  HHlevelDat <- CounterFactualResultsPost$ExplodedData[, .(
    p_1 = ifelse(mode1(alpha - psi) > 0, 1, 0),
    p_2 = ifelse(mode1(alpha - psi) > alphampsicutoff, 1, 0),
    psi = mode1(psi),
    alpha = mode1(alpha),
    omega = mode1(omega),
    subs = mode1(subs),
    best_guess_AV_old = mode1(best_guess_AV_old) / 100,
    best_guess_sumrate = mode1(best_guess_sumrate) / 100,
    best_guess_netprem_old = mode1(best_guess_netprem_old) / 100,
    grossprem_old = mode1(grossprem_old) / 100,
    ratingfactor = mode1(ratingfactor),
    InPreSample = mode1(!is.na(grossprem_old))),
    by = c("hhid", "year")]
  HHlevelDat[, Market := ifelse(hhid < 0, "SmallGroup", "Individual")]
  
  SGpretab <- with(HHlevelDat[Market == "SmallGroup" & InPreSample == T], c(
    mean(pmin(1 / alpha, costwinsorization), na.rm = T),
    mean(pmin((best_guess_AV_old + omega * (best_guess_AV_old ^ 2)) / alpha, costwinsorization) - pmin(best_guess_AV_old / alpha, costwinsorization), na.rm = T),
    0.000, # share uninsured
    sum(best_guess_AV_old == 0.6, na.rm = T) / sum(!is.na(best_guess_AV_old)),
    sum(best_guess_AV_old == 0.7, na.rm = T) / sum(!is.na(best_guess_AV_old)),
    sum(best_guess_AV_old == 0.8, na.rm = T) / sum(!is.na(best_guess_AV_old)),
    sum(grossprem_old / ratingfactor * as.numeric(best_guess_AV_old == 0.6), na.rm = T) / sum(best_guess_AV_old == 0.6, na.rm = T),
    sum(grossprem_old / ratingfactor * as.numeric(best_guess_AV_old == 0.7), na.rm = T) / sum(best_guess_AV_old == 0.7, na.rm = T),
    sum(grossprem_old / ratingfactor * as.numeric(best_guess_AV_old == 0.8), na.rm = T) / sum(best_guess_AV_old == 0.8, na.rm = T),
    sum((grossprem_old * (1 - nu) * (1 - best_guess_sumrate)  / ratingfactor) * as.numeric(best_guess_AV_old == 0.6), na.rm = T) / sum(best_guess_AV_old == 0.6, na.rm = T),
    sum((grossprem_old * (1 - nu) * (1 - best_guess_sumrate) / ratingfactor) * as.numeric(best_guess_AV_old == 0.7), na.rm = T) / sum(best_guess_AV_old == 0.7, na.rm = T),
    sum((grossprem_old * (1 - nu) * (1 - best_guess_sumrate) / ratingfactor) * as.numeric(best_guess_AV_old == 0.8), na.rm = T) / sum(best_guess_AV_old == 0.8, na.rm = T),
    mean(best_guess_sumrate * grossprem_old, na.rm = T), #gvmt spend
    ifelse(subsidytypeSG == "Ratio", 
      mean(nu * (1 - best_guess_sumrate) * grossprem_old, na.rm = T), 
      mean(subs)), # empl spend
    sum(!is.na(grossprem_old))
  ))
  
  # Filling in table:
  print("Making Table...")
  Table1$IndPre <- c(
    lapply(pretab[Market == "Individual", .(CS_1, p_1, CS_2, p_2, ExpLambda, ExpMHS, PUninsured, PBronze, PSilver, PGold, AveBronzeBaseGross, AveSilverBaseGross, AveGoldBaseGross, AveBronzePremNet, AveSilverPremNet, AveGoldPremNet, GovmntSpend), ], median, na.rm = T),
    NA, nrow(pretab[Market == "Individual"]))
  Table1$IndPost <- c(
    lapply(posttab[Market == "Individual", .(CS_1, p_1, CS_2, p_2, ExpLambda, ExpMHS, PUninsured, PBronze, PSilver, PGold, AveBronzeBaseGross, AveSilverBaseGross, AveGoldBaseGross, AveBronzePremNet, AveSilverPremNet, AveGoldPremNet, GovmntSpend), ], median, na.rm = T),
    NA, nrow(posttab[Market == "Individual"]))
  try(Table1$IndChange <- as.numeric(Table1$IndPost) - as.numeric(Table1$IndPre))
  Table1$SGPre <- c(NA, NA, NA, NA, SGpretab)
  Table1$SGPreSamplePost <- c(
    lapply(posttab[Market == "SmallGroup" & InPreSample == T, .(CS_1, p_1, CS_2, p_2, ExpLambda, ExpMHS, PUninsured, PBronze, PSilver, PGold, AveBronzeBaseGross, AveSilverBaseGross, AveGoldBaseGross, AveBronzePremNet, AveSilverPremNet, AveGoldPremNet, GovmntSpend, EmployerSpend), ], median, na.rm = T),
    sum(rev(SGpretab)[1]))
  try(Table1$SGChange <- as.numeric(Table1$SGPreSamplePost) - as.numeric(Table1$SGPre))
  Table1$SGPost <- c(
    lapply(posttab[Market == "SmallGroup", .(CS_1, p_1, CS_2, p_2, ExpLambda, ExpMHS, PUninsured, PBronze, PSilver, PGold, AveBronzeBaseGross, AveSilverBaseGross, AveGoldBaseGross, AveBronzePremNet, AveSilverPremNet, AveGoldPremNet, GovmntSpend, EmployerSpend), ], median, na.rm = T),
    nrow(posttab[Market == "SmallGroup"]))
  
  print(paste0(
    "In the small group pre-sample, ",
    sum(HHlevelDat[Market == "SmallGroup"]$best_guess_AV_old == 0.9, na.rm = T),
    " households buy a platinum plan."
  ))
  print(paste0(
    "In the individual market pre-sample, there are ",
    mean(CounterFactualResultsPre$ExplodedData[hhid > 0, .(nchoice = .N), by = c("hhid", "year")]$nchoice, na.rm = T),
    " average number of plans in choice set."
  ))
  
  return(Table1)
}


flipdim <- function(mat){
  if(nrow(mat) != ncol(mat)){ stop("Must be square") }
  matnew <- matrix(nrow = nrow(mat), ncol = ncol(mat))
  for(i in 1:nrow(mat)){
    for(j in 1:ncol(mat)){
      matnew[i,j] <- mat[j,i]
    }
  }
  rownames(matnew) <- rownames(mat)
  colnames(matnew) <- colnames(mat)
  return(matnew)
}


CreateTransitionTables <- function(
  CounterFactualResultsPre,
  CounterFactualResultsPost,
  costwinsorization = 89.37,
  SGpremsubs = F, 
  subsidytype = "Ratio"
){
  
  # Labeling Transition type:
  CounterFactualResultsPre$ExplodedData[, 
    type := fcase(x_av == 0, "0Uninsured",
      x_av == 0.6, "6Bronze",
      x_av == 0.8, "8Gold", 
      default = "7Silver"), ]
      
  CounterFactualResultsPost$ExplodedData[, 
    type := fcase(x_av == 0, "0Uninsured",
      x_av == 0.6, "6Bronze",
      x_av == 0.8, "8Gold", 
      default = "7Silver"), ]
  
  # Collapsing to type:
  CounterCollapsedPre <- CounterFactualResultsPre$ExplodedData[, .(
    alpha = mean(alpha), 
    share = sum(share)), 
    by = c("year", "hhid", "type")]
    
  CounterCollapsedPost <- CounterFactualResultsPost$ExplodedData[, .(
    alpha = mean(alpha), 
    share = sum(share)), 
    by = c("year", "hhid", "type")]
  
  # Merging Pre and Post Tables (allowing for all possible type combinations): 
  CounterFactualResultsMerged <- merge(
    CounterCollapsedPre, 
    CounterCollapsedPost, 
    by = c("year", "hhid"), 
    all.x = T, 
    suffixes = c(".pre", ".post"), 
    allow.cartesian = T)
    
  # Create Transition Tables:
  CounterFactualResultsMerged$n_withshare = sum(CounterFactualResultsMerged$share.pre * CounterFactualResultsMerged$share.post)
  
  TransitionDT <- CounterFactualResultsMerged[, .(
    ExpLambda = sum(share.pre * share.post * pmin(1 / alpha.pre, costwinsorization)) / sum(share.pre * share.post),
    P_transition = sum(share.pre * share.post / n_withshare)),
    by = c("type.pre", "type.post")]
    
  P_transition_table <- dcast(TransitionDT, formula = type.pre ~ type.post, value.var = "P_transition")
  P_transition_mat <- as.matrix(P_transition_table[, 2:5])
  rownames(P_transition_mat) <- P_transition_table$type.pre
  ExpLambda_transition_table <- dcast(TransitionDT, formula = type.pre ~ type.post, value.var = "ExpLambda")
  ExpLambda_transition_mat <- as.matrix(ExpLambda_transition_table[, 2:5])
  rownames(ExpLambda_transition_mat) <- ExpLambda_transition_table$type.pre
  
  # Differences Tables:
  P_transition_diff <- (P_transition_mat - flipdim(P_transition_mat))
  ExpLambda_transition_diff <- (ExpLambda_transition_mat - flipdim(ExpLambda_transition_mat))
  
  return(list(probs_mat = P_transition_mat, 
    spend_mat = ExpLambda_transition_mat, 
    probs_diff = P_transition_diff, 
    spend_diff = ExpLambda_transition_diff))
}

generateOutcomeColumns <- function(
  CounterFactualResults,
  costwinsorization = 89.37,
  SGpremsubs = F, # Are SG households given g'ment subsidies?
  subsidytypeSG = "Ratio", # Other options are "Fixed" or "CSCAdjust"
  subsidytypeInd = "Ratio" # Other options are "Fixed" or "CSCAdjust"
){

  print("Creating HH level results.")
  posttab <- CounterFactualResults$ExplodedData[, .(
    acg = mode1(sum_concurrent_risk),
    CS = mode1(1 / (alpha - psi)) * log(sum(exp(V))),
    p_1 = ifelse(mode1(alpha - psi) > 0, 1, 0),
    p_2 = ifelse(mode1(alpha - psi) > alphampsicutoff, 1, 0),
    PUninsured = sum(share * as.numeric(x_av == 0)),
    PBronze = sum(share * as.numeric(x_av == 0.6)), 
    PSilver = sum(share * as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94))),
    PGold = sum(share * as.numeric(x_av == 0.8)),
    AveBronzePremNet = sum(as.numeric(x_av == 0.6) * ifelse(is.na(price_hh), 0, price_hh) * share) / sum(as.numeric(x_av == 0.6) * share), 
    AveSilverPremNet = sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94)) * ifelse(is.na(price_hh), 0, price_hh) * share) / sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94)) * share),
    AveGoldPremNet = sum(as.numeric(x_av == 0.8) * ifelse(is.na(price_hh), 0, price_hh) * share) / sum(as.numeric(x_av == 0.8) * share), 
    AveBronzeSubsidy = sum(as.numeric(x_av == 0.6) * (ratingfactor * price - ifelse(is.na(price_hh), 0, price_hh)) * share) / sum(as.numeric(x_av == 0.6) * share), 
    AveSilverSubsidy = sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94)) * (ratingfactor * price - ifelse(is.na(price_hh), 0, price_hh)) * share) / sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94)) * share), 
    AveGoldSubsidy = sum(as.numeric(x_av == 0.8) * (ratingfactor * price - ifelse(is.na(price_hh), 0, price_hh)) * share) / sum(as.numeric(x_av == 0.8) * share), 
    AveBronzeBaseGross = sum(as.numeric(x_av == 0.6) * ifelse(is.na(price), 0, price) * share) / sum(as.numeric(x_av == 0.6) * share), 
    AveSilverBaseGross = sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94)) * ifelse(is.na(price), 0, price) * share) / sum(as.numeric(x_av %in% c(0.7, 0.73, 0.87, 0.94)) * share), 
    AveGoldBaseGross = sum(as.numeric(x_av == 0.8) * ifelse(is.na(price), 0, price) * share) / sum(as.numeric(x_av == 0.8) * share), 
    GovmntSpend = ifelse(SGpremsubs == T | hhid > 0, 
      ifelse(subsidytypeInd == "Ratio",
        sum(share * subs * ratingfactor * price), 
        sum(share * ifelse(subs < ratingfactor * price, subs, ratingfactor * price))),
      sum(share * cmax((best_guess_sumrate / 100) * ratingfactor * price, nu * (best_guess_sumrate  / 100) * (grossprem_old / 100)))),
    EmployerSpend = ifelse(SGpremsubs == T | hhid > 0, 
      0,
      ifelse(subsidytypeSG == "Ratio", 
        sum(share * (share / ( 1 - share)) * price_hh), 
        subs)),
    InPreSample = mode1(!is.na(grossprem_old))
  ),
  by = c("hhid", "year")]
  
  posttab[, Market := ifelse(hhid < 0, "SmallGroup", "Individual")]
  posttab[, CS_1 := ifelse(p_1 == 1, CS, NA), ]
  posttab[, CS_2 := ifelse(p_2 == 1, CS, NA), ]
  posttab[, CS_2_highacg := ifelse(p_2 == 1 & acg > acgHL_boundary, CS, NA), ]
  posttab[, CS_2_lowacg := ifelse(p_2 == 1 & acg <= acgHL_boundary, CS, NA), ]
  
  print(paste0("Ind Mkt N: ", nrow(posttab[Market == "Individual"])))
  print(paste0("SG Mkt N: ", nrow(posttab[Market == "SmallGroup"])))

  print("Creating HH level results.")
  Results_List <- list()
  
  # Set variables to use
  welfare_vars_list <- list(ind_mkt = c("CS_1", "CS_2", "GovmntSpend", "CS_2_lowacg", "CS_2_highacg"), 
    sg_mkt = c("CS_1", "CS_2", "GovmntSpend", "EmployerSpend", "CS_2_lowacg", "CS_2_highacg"))
  mktshr_vars <- c("PUninsured", "PBronze", "PSilver", "PGold")
  grossprem_vars <- c("AveBronzeBaseGross", "AveSilverBaseGross", "AveGoldBaseGross")
  netprem_vars <- c("AveBronzePremNet", "AveSilverPremNet", "AveGoldPremNet") 
  subsidy_vars <- c("AveBronzeSubsidy", "AveSilverSubsidy", "AveGoldSubsidy")
  
  markets <- list("Individual", "SmallGroup", c("Individual", "SmallGroup"))
  for(i in 1:2){
    Results_List[[i]] <- list()
    Results_List[[i]]$welfare <- unlist(lapply(posttab[Market %in% markets[[i]], welfare_vars_list[[i]], with = F], mean, na.rm = T))
    Results_List[[i]]$mktshr <- unlist(lapply(posttab[Market %in% markets[[i]], mktshr_vars, with = F], mean, na.rm = T))
    Results_List[[i]]$grossprem <- with(posttab[Market %in% markets[[i]]], 
      c(weighted.mean(AveBronzeBaseGross, PBronze, na.rm = T), 
      weighted.mean(AveSilverBaseGross, PSilver, na.rm = T),
      weighted.mean(AveGoldBaseGross, PGold, na.rm = T)))
    Results_List[[i]]$netprem <- with(posttab[Market %in% markets[[i]]], 
      c(weighted.mean(AveBronzePremNet, PBronze, na.rm = T), 
      weighted.mean(AveSilverPremNet, PSilver, na.rm = T),
      weighted.mean(AveGoldPremNet, PGold, na.rm = T)))
    Results_List[[i]]$subsidy <- with(posttab[Market %in% markets[[i]]], 
      c(weighted.mean(AveBronzeSubsidy, PBronze, na.rm = T), 
      weighted.mean(AveSilverSubsidy, PSilver, na.rm = T),
      weighted.mean(AveGoldSubsidy, PGold, na.rm = T)))
    }

  return(Results_List)
}


generateSGOutcomeColumnsPre <- function(
  CounterFactualResults,
  costwinsorization = 89.37,
  SGpremsubs = F, 
  subsidytypeSG = "Ratio", # Other options are "Fixed" or "CSCAdjust"
  subsidytypeInd = "Ratio" # Other options are "Fixed" or "CSCAdjust"
){

  print("Creating HH level results.")
  HHlevelDat <- CounterFactualResults$ExplodedData[, .(
    p_1 = ifelse(mode1(alpha - psi) > 0, 1, 0),
    p_2 = ifelse(mode1(alpha - psi) > alphampsicutoff, 1, 0),
    acg = mode1(sum_concurrent_risk),
    psi = mode1(psi),
    alpha = mode1(alpha),
    omega = mode1(omega),
    subs = mode1(subs),
    best_guess_AV_old = mode1(best_guess_AV_old) / 100,
    best_guess_sumrate = mode1(best_guess_sumrate) / 100,
    best_guess_netprem_old = mode1(best_guess_netprem_old) / 100,
    grossprem_old = mode1(grossprem_old) / 100,
    ratingfactor = mode1(ratingfactor),
    InPreSample = mode1(!is.na(grossprem_old))),
    by = c("hhid", "year")]
  HHlevelDat[, Market := ifelse(hhid < 0, "SmallGroup", "Individual")]
  
  foo <- HHlevelDat[Market == "SmallGroup" & InPreSample == T]
  print(paste0("SG Base N:", nrow(foo)))

  print("Creating HH level results.")
  Results_List <- list()
  
  # Set variables to use
  welfare_vars_list <- list(ind_mkt = c("CS_1", "CS_2", "GovmntSpend", "CS_2_lowacg", "CS_2_highacg"), 
    sg_mkt = c("CS_1", "CS_2", "GovmntSpend", "EmployerSpend", "CS_2_lowacg", "CS_2_highacg"))
  mktshr_vars <- c("PUninsured", "PBronze", "PSilver", "PGold")
  grossprem_vars <- c("AveBronzeBaseGross", "AveSilverBaseGross", "AveGoldBaseGross")
  netprem_vars <- c("AveBronzePremNet", "AveSilverPremNet", "AveGoldPremNet") 
  subsidy_vars <- c("AveBronzeSubsidy", "AveSilverSubsidy", "AveGoldSubsidy")
    
  Results <- list()

  Results$welfare <- with(foo, 
    c(
      NA, 
      NA, 
      mean(best_guess_sumrate * grossprem_old, na.rm = T), 
      mean(nu * (1 - (best_guess_sumrate)) * grossprem_old, na.rm = T), 
      NA, 
      NA))

  Results$mktshr <- with(foo, 
    c(0.000, # share uninsured
    sum(best_guess_AV_old == 0.6, na.rm = T) / sum(!is.na(best_guess_AV_old)),
    sum(best_guess_AV_old %in% c(0.7, 0.73, 0.87), na.rm = T) / sum(!is.na(best_guess_AV_old)),
    sum(best_guess_AV_old == 0.8, na.rm = T) / sum(!is.na(best_guess_AV_old))))

  Results$grossprem <- with(foo, 
    c(sum(grossprem_old / ratingfactor * as.numeric(best_guess_AV_old == 0.6), na.rm = T) / sum (best_guess_AV_old == 0.6, na.rm = T),
    sum(grossprem_old / ratingfactor * as.numeric(best_guess_AV_old %in% c(0.7, 0.73, 0.87)), na.rm = T) / sum(best_guess_AV_old %in% c(0.7, 0.73, 0.87), na.rm = T),
    sum(grossprem_old / ratingfactor * as.numeric(best_guess_AV_old == 0.8), na.rm = T) / sum(best_guess_AV_old == 0.8, na.rm = T)))

  Results$netprem <- with(foo, 
    c(sum((grossprem_old * (1 - nu) * (1 - best_guess_sumrate)) * as.numeric(best_guess_AV_old == 0.6), na.rm = T) / sum(best_guess_AV_old == 0.6, na.rm = T),
    sum((grossprem_old * (1 - nu) * (1 - best_guess_sumrate)) * as.numeric(best_guess_AV_old %in% c(0.7, 0.73, 0.87)), na.rm = T) / sum(best_guess_AV_old %in% c(0.7, 0.73, 0.87), na.rm = T),
    sum((grossprem_old * (1 - nu) * (1 - best_guess_sumrate)) * as.numeric(best_guess_AV_old == 0.8), na.rm = T) / sum(best_guess_AV_old == 0.8, na.rm = T)))

  Results$subsidy <- with(foo, 
    c(sum((grossprem_old * (1 - ((1 - nu) * (1 - best_guess_sumrate)))) * as.numeric(best_guess_AV_old == 0.6), na.rm = T) / sum (best_guess_AV_old == 0.6, na.rm = T),
    sum((grossprem_old * (1 - ((1 - nu) * (1 - best_guess_sumrate)))) * as.numeric(best_guess_AV_old %in% c(0.7, 0.73, 0.87)), na.rm = T) / sum(best_guess_AV_old %in% c(0.7, 0.73, 0.87), na.rm = T),
    sum((grossprem_old * (1 - ((1 - nu) * (1 - best_guess_sumrate)))) * as.numeric(best_guess_AV_old == 0.8), na.rm = T) / sum(best_guess_AV_old == 0.8, na.rm = T)))

  return(Results)
}
